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Abstract 

We introduce a level set based approach to Bayesian geometric inverse problems. In these 
problems the interface between different domains is the key unknown, and is realized as the 
level set of a function. This function itself becomes the object of the inference. Whilst the 
level set methodology has been widely used for the solution of geometric inverse problems, the 
Bayesian formulation that we develop here contains two significant advances: firstly it leads to 
a well-posed inverse problem in which the posterior distribution is Lipschitz with respect to the 
observed data; and secondly it leads to computationally expedient algorithms in which the level 
set itself is updated implicitly via the MCMC methodology applied to the level set function - 
no explicit velocity field is required for the level set interface. Applications are numerous and 
include medical imaging, modelling of subsurface formations and the inverse source problem; 
our theory is illustrated with computational results involving the last two applications. 


1 Introduction 

Geometric inverse problems, in which the interfaces between different domains are the primary un¬ 
known quantities, are ubiquitous in applications including medical imaging problems such as EIT 
[9] and subsurface flow [5]; they also have an intrinsically interesting mathematical structure [30], 
In many such applications the data is sparse, so that the problem is severely under-determined, and 
noisy. For both these reasons the Bayesian approach to the inverse problem is attractive as the proba¬ 
bilistic formulation allows for regularization of the under-determined, and often ill-posed, inversion 
via the introduction of priors, and simultaneously deals with the presence of noise in the observa¬ 
tions [32, 45]. The level set method has been a highly successful methodology for the solution of 
classical, non-statistical, inverse problems for interfaces since the seminal paper of Santosa [43]; 
see for example [15, 11, 12, 47, 33, 34, 46, 13, 18, 19, 48, 28, 3, 6] and for related Bayesian level 
set approaches see [49, 36, 37, 40]. 
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In this paper we marry the level set approach with the Bayesian approach to geometric inverse 
problems. This leads to two significant advances: firstly it leads to a well-posed inverse problem 
in which the posterior distribution is Lipschitz with respect to the observed data, in the Hellinger 
metric - there is no analogous well-posedness theory for classical level set inversion; and secondly 
it leads to computationally expedient algorithms in which the interfaces are updated implicitly via 
the Markov chain Monte Carlo (MCMC) methodology applied to the level set function - no ex¬ 
plicit velocity field is required for the level set interface. We highlight that the recent paper [48] 
demonstrates the potential for working with regularized data misfit minimization in terms of a level 
set function, but is non-Bayesian in its treatment of the problem, using instead simulated annealing 
within an optimization framework. On the other hand the paper [49] adopts a Bayesian approach 
and employs the level set method, but requires a velocity field for propagation of the interface and 
does not have the conceptual and implementational simplicity, as well as the underlying theoretical 
basis, of the method introduced here. The papers [ 36 , 37 , 40], whilst motivated by the Bayesian 
approach, use the ensemble Kalman filter and are therefore not strictly Bayesian - the method does 
not deliver a provably reliable approximation of the posterior distribution except for linear Gaussian 
inverse problems. 

The key idea which underpins our work is this. Both the theory and computational practice of 
the level set method for geometric inverse problems is potentially hampered by the fact that the 
mapping from the space of the level set function to the physical parameter space is discontinuous. 
This discontinuity occurs when the level set function is flat at the critical levels, and in particular 
where the desired level set has non-zero Lebesgue measure. This is dealt with in various ad hoc 
ways in the applied literature. The beauty of the Bayesian approach is that, with the right choice 
of prior in level set space, these discontinuities have probability zero. As a result a well-posedness 
theory (the posterior is Lipschitz in the data) follows automatically, and computational algorithms 
such as MCMC may be formulated in level set space. We thus have practical algorithms which are 
simultaneously founded on a sound theoretical bedrock. 

In section 2 we set up the inverse problem of interest, describe the level set formulation, and 
state assumptions under which we have a well-posed inverse problem for the level set function. 
We also characterize the discontinuity set of the level set map and demonstrate the existence of 
Gaussian priors for which this discontinuity set is a probability zero event. In section 3 we describe 
two examples - inverse gravimetry and permeability determination in groundwater flow - which can 
be shown to satisfy the theoretical framework of the preceding section and hence for which there is 
a well-posed inverse problem for the level set function. Section 4 contains numerical experiments 
for both of these examples, demonstrating the potential of the methodology, and also highlighting 
important questions for future research. We conclude in section 5, and then the two appendices 
contain various technical details and proofs which have been deferred in order to maintain the flow 
of ideas in the main body of the article. 
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2 Bayesian Level Set Inversion 

2.1 The Inverse Problem 

This paper is concerned with inverse problems of the following type: recover function k G X := 
L q (D\ M), D a bounded open set in M 2 , from a finite set of noisily observed linear functionals 
{Oj }j =1 of p £ V, for some Banach space V, where p = G(u) for nonlinear operator G G C(X , V ) . 
Typically, for us, k will represent input data for a partial differential equation (PDE), p the solution 
of the PDE and G the solution operator mapping the input n to the solution p. Collecting the linear 
functionals into a single operator O : V —> M J and assuming additive noise r\ we obtain the inverse 
problem of finding n from y where 

y = (O oG)(k) +r). ( 2 . 1 ) 

This under-determined inverse problem is well-adapted to both the classical [ 20 ] and Bayesian [ 17 ] 
approaches to regularized inversion, because O o G e C(X. M J ). 

2.2 Level Set Parameterization 

There are many inverse problems where k is known a priori to have the form 

n 

n( x ) = ^2nil Di (x)] (2.2) 

1=1 

here Id denotes the indicator function of subset B c M 2 , {A:}” = i are subsets of D such that 
U" =1 _Dj = D and D, n D 3 = 0, and the {«:«}”= i are known positive constants. 1 In this setting 
the Di become the primary unknowns and the level set method is natural. Given integer n fix the 
constants q £ R. for i = 0, • • • , n with — oo = c 0 < c\ < • • • < c n = oo and consider a real-valued 
continuous level set function u : D —y R. We can then define the £>* by 

Di = {x e D \ q_i < u{x) < q}. (2.3) 

It follows that Di fi Dj = 0 for i,j > 1, i f j. For later use define the z-th level set D ( - = 
Di fl D l+ i = {x G D u{x) = q). Let U = C(D; M) and, given the positive constants {/q}” =1 , we 
define the level set map F : U —> X by 

n 

(■ Fu ) (x) -> «:(a:) = ^2 K i (x). (2.4) 

2=1 

We may now define Q = O o G o F :[/—)■ M J and reformulate the inverse problem in terms of the 
level set function u: find u from y where 

V = Q{u) + p. (2.5) 

'Generalization to the Ki being unknown constants, or unknown smooth functions on each domain Di, are possible 
but will not be considered explicitly in this paper. Our focus is on the geometry of the interfaces implied by the Dj. 
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However, because F : U —» X, and hence Q : U —» R J , is discontinuous, the classical regular¬ 
ization theory for this form of inverse problem is problematic. Whilst the current state of the art 
for Bayesian regularization assumes continuity of Q for inverse problems of the form (2.5), we will 
demonstrate that the Bayesian setting can be generalized to level set inversion. This will be achieved 
by a careful understanding of the discontinuity set for F, and an understanding of probability mea¬ 
sures for which this set is a measure zero event. 


2.3 Well-Posed Bayesian Inverse Problem 

In the Bayesian approach to inverse problems, all quantities in (2.5) are treated as random variables. 
Let U denote a separable Banach space and define a complete 2 probability space ( U , E, // 0 ) for the 
unknown u. (In our applications U will be the space C(D: R) but we state our main theorem in 
more generality). Assume that the noise 77 is a random draw from the centered Gaussian Q 0 := 
A/"(0, T); we also assume that 77 and u are independent . 3 We may now define the joint random 
variable (u, y) £ U x R J . The Bayesian approach to the inverse problem (2.5) consists of seeking 
the posterior probability distribution y y of the random variable u\y, u given y. This measure jt v 
describes our probabilistic knowledge about u on the basis of the measurements y given by (2.5) 
and the prior information /i 0 on u. 

We define the least squares function 4> : U x R J —y R + by 

$( u ',y) = ^\y~S(u)\l ( 2 . 6 ) 

with | • | := |r - ^ . |. We now state a set of assumptions under which the posterior distribution is 
well-defined via its density with respect to the prior distribution, and is Lipschitz in the Hellinger 
metric, with respect to data y. 

Assumptions 2.1. The least squares function <E> : U x R J —» R and probability measure // 0 on the 
measure space ( U , E) satisfy the following properties: 

1. for every r > 0 there is a K = K(r) such that, for all a £ U and all y <G M J with |y| r < r, 

0 < $(m; y) < K] 

2 . for any fixed y £ R J , $(-;y) : U —* R, is continuous 7 i 0 -almost surely on the complete 
probability space ([/, E, y 0 ); 

3. for yi, 7/2 £ R J with max{|y!|r, 1 2 / 2 1r} < f, there exists a C = C(r) such that, for all u £ U, 

\$(u-,yi) - $(u;y 2 )\ < C\yi - y 2 \r- 

2 Complete probability space is defined at the start of Appendix 2 

3 Allowing for non-Gaussian 77 is also possible, as is dependence between 77 and u; however we elect to keep the 
presentation simple. 
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The Hellinger distance between p and /// is defined as 


^Hell(/L AO 
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for any measure v with respect to which p and p' are absolutely continuous; the Hellinger distance 
is, however, independent of which reference measure v is chosen. We have the following: 

Theorem 2.2. Assume that the least squares function $ : U x R J —>■ E and the probability measure 
Po on the measure space (U. E) satisfy Assumptions 2.1. Then p y IM> with Radon-Nikodym 
derivative 

= tf exp (—$(«; y)) (2.7) 

whera, for y almost surely, 

Z := exp(-$(u;y))p 0 (d,u) > 0 . 

Ju 

Furthermore p y is locally Lipschitz with respect to y, in the Hellinger distance: for all y, y' with 

max{|y| r , 12/' | r} < r, there exists aC = C(r) > 0 such that 

dne\\{p y ,P V ) < C\y — y' |r- 


This implies that, for all / e L‘f i0 (U : S) for separable Banach space S', 


\\W v f{u)-W vl f(u)\\s<C\y-y' 


Remarks 2.3. • The interpretation of this result is very natural, linking the Bayesian picture 

with least squares minimisation: the posterior measure is large on sets where the least squares 
function is small, and vice-versa, all measured relative to the prior p Q . 


• The key technical advance in this theorem over existing theories overviewed in [17] is that 
$(•;?/) is only continuous p 0 —almost surely; existing theories typically use that <!>(•; y) is 
continuous everywhere on U and that p 0 (U) = 1; these existing theories cannot be used 
in the level set inverse problem, because of discontinuities in the level set map. Once the 
technical Lemma 6.1 has been established, which uses p 0 —almost sure continuity to establish 
measurability, the proof of the theorem is a straightforward application of existing theory; we 
therefore defer it to Appendix 1. 


• What needs to be done to apply this theorem in our level set context is to identify the sets 
of discontinuity for the map Q, and hence $(•: ;y), and then construct prior measures p 0 for 
which these sets have measure zero. We study these questions in general terms in the next two 
subsections, and then, in the next section, demonstrate two test model PDE inverse problems 
where the general theory applies. 
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• The consequences of this result are wide-ranging, and we name the two primary ones: firstly 
we may apply the mesh-independent MCMC methods overviewed in [16] to sample the poste¬ 
rior distribution efficiently; and secondly the well-posedness gives desirable robustness which 
may be used to estimate the effect of other perturbations, such as approximating G by a nu¬ 
merical method, on the posterior distribution [17]. 

2.4 Discontinuity Sets of F 

We return to the specific setting of the level set inverse problem where U = C(D; E). We first note 
that the level set map F : U —> L°°(D ; E) is not continuous except at points u E U where no level 
crossings occur in an open neighbourhood. However as a mapping F : U —>■ L q (D ; E) for q < oc 
the situation is much better: 

Proposition 2.4. For u E C{D) and 1 < q < oo, the level set map F : C(D) —> L q (D) is 
continuous at u if and only ifm(Df) = 0 for all i = 1, • • • ,7i — l. 

The proof is given in Appendix 1. For the inverse gravimetry problem considered in the next 
section the space X is naturally L 2 (D ; E) and we will be able to directly use the preceding proposi¬ 
tion to establish almost sure continuity of F and hence Q. For the groundwater flow inverse problem 
the space X is naturally L°°(D ; E) and we will not be able to use the proposition in this space to 
establish almost sure continuity of F. However, we employ recent Lipschitz continuity results for 
G on L q (D ; E), q < oc to establish the almost sure continuity of Q [8]. 

2.5 Prior Gaussian Measures 

Let D denote a bounded open subset of E 2 . For our applications we will use the following two con¬ 
structions of Gaussian prior measures go which are Gaussian A^(0, C, ), i = 1, 2 on Hilbert function 
space Hi, i = 1,2. 

• Af(0, Ci) on 

Hi := {u : D —y E | u £ L 2 (D ; E)), / u(x ) dx = 0}, 

Jd 

where 

Ci = A~ a and A := -A (2.8) 

with domain 4 

D(A) := {u : D —> E | u E H 2 (D ; E), Xu ■ v = 0 on dD and / u(x) dx = 0}. 

Jd 

4 Here v denotes the outward normal. 
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• JV(0, C 2 ) on H 2 := L 2 (D ; M) with C 2 : 77 2 —> H being the integral operator 

C 2 (j)(x) = J c(x, y)<f>(y) dy with c(x, y ) = exp ^ (2.9) 

In fact, in the inverse model arising from groundwater flow studied in [26, 25], the Gaussian measure 
A/"(0, C\) was taken as the prior measure for the logarithm of the permeability. On the other hand the 
Gaussian measure 7V*(0, C 2 ) is widely used to model the earth’s subsurface [39] as draws from this 
measure generate smooth random functions in which the parameter L sets the spatial correlation 
length. For both of these measures it is known that, under suitable conditions on the domain D, 
draws are almost surely in C(D\ M); see [17], Theorems 2.16 and 2.18 for details. 

Assume that a > 1 in (2.8), then the Gaussian random function with measure // 0 defined in 
either case above has the property that, for U := C'(D\ R), /i 0 ( U) = 1. Since U is a separable 
Banach space /v, 0 can be redefined as a Gaussian measure on U. Furthermore it is possible to define 
the appropriate cr-algebra E in such a way that (U, E, /.i 0 ) is a complete probability space; for details 
see Appendix 2. We have the following, which is a subset of what is proved in Proposition 7.2. 

Proposition 2.5. Consider a random function u drawn from one of the Gaussian probability mea¬ 
sures fj, 0 on U given above. Then rri(Df) = 0, /.(o-almost surely, for % = 1. • • • , n. 

This, combined with Proposition 2.4, is the key to making a rigorous well-posed formulation of 
Bayesian level set inversion. Together the results show that priors may be constructed for which the 
problematic discontinuities in the level set map are probability zero events. In the next section we 
demonstrate how the theory may be applied, by consider two examples. 

3 Examples 

3.1 Test Model 1 (Inverse Potential Problem) 

Let D C M 2 be a bounded open set with Lipschitz boundary. Consider the PDE 

A p = K in D, p = 0 on dD. (3.1) 

If k G X := L 2 (D) it follows that there is a unique solution p £ Hq(D). Furthermore A p £ L 2 (D), 
so that the Neumann trace can be defined in V := H~^(dD) by the following Green’s formula: 

= [ Ap<pdx+ [ VpVpdx 
\dv / dD J D J D 

for ip £ H l (D). Here v is the unit outward normal vector on dD and (-,-)dD denotes the dual 
pairing on the boundary. We can then define the bounded linear map G : X —» V by G(k) = 4^. 
Now assume that the source term k has the form 


k(x) = I dA x ) 
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for some D\ C D. The inverse potential problem is to reconstruct the support D\ from measure¬ 
ments of the Neumann data of p on OD. In the case where the Neumann data is measured every¬ 
where on the boundary dD, and where the domain D y is assumed to be star-shaped with respect to 
its center of gravity, the inverse problem has a unique solution; see [30, 31] for details of this theory 
and see [23, 13] for discussion of numerical methods for this inverse problem. We will study the 
underdetermined case where a finite set of bounded linear functionals Oj : V —> M are measured, 
noisily, on dD: 



Concatenating we have y = (O o G)(k) + //. Representing the boundary of D\ as the zero level set 
of function u & U := C(D: M) we write the inverse problem in the form (2.5): 

y — (O o G o F)(u) + rj. (3.3) 

Since multiplicity and uncertainty of solutions are then natural, we will adopt a Bayesian approach. 

Notice that the level set map F : U — >■ X is bounded: for all u G U we have ||F(u) || y < 
Vol(-D) := f D dx. Since G : X —» V and O : V —y M J are bounded linear maps it follows that 
Q = OoGoF:U R J is bounded: we have constant C + G M + such that, for all u G U, 
\G(u) | < C + . From this fact Assumptions 2.1(i) and (iii) follow automatically. Since both G : X —>■ 

V and O : V —» M J are bounded, and hence continuous, linear maps, the discontinuity set of Q 
is determined by the the discontinuity set of F : U —» X. By Proposition 2.4 this is precisely the 
set of functions for which the measure of the level set { u(x ) = 0} is zero. By Proposition 2.5 this 
occurs with probability zero for both of the Gaussian priors specified there and hence Assumptions 
2.1(ii) holds with these priors. Thus Theorem 2.2 applies and we have a well-posed Bayesian inverse 
problem for the level set function. 

3.2 Test Model 2 (Discontinuity Detection in Groundwater Flow) 

Consider the single-phase Darcy-flow model given by 

— V • (kVp) = / inD,p = 0 on dD. (3.4) 

Here D is a bounded Lipschitz domain in R 2 , k the real-valued isotropic permeability function and 
p the fluid pressure. The right hand side / accounts for the source of groundwater recharge. Let 

V = H'(D; M), X = L°°(D ; M) and V* denote the dual space of F. If / G V* and A+ := {k G 
X : essinf^gD/^x) > ft min > 0} then G : X + V defined by G(k) = p is Lipschitz continuous 
and 

||G r («)||v = |b|k ^ ||/||v*/^mm- (3.5) 

We consider the practically useful situation in which the permeability function k is modelled as 
piecewise constant on different regions {A}[* =1 whose union comprise D; this is a natural way to 
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characterize heterogeneous subsurface structure in a physically meaningful way. We thus have 


n 

k(x) = y^KiIp.(x) 

2=1 

where {A}™= 1 are subsets of D such that U'' =1 A = D and A n Dj = 0, and where the {ftj}" =1 
are positive constants. We let K rnitl = min, k, . 

Unique reconstruction of the permeability in some situations is possible if the pressure p is mea¬ 
sured everywhere [2, 41]. The inverse problem of interest to us is to locate the discontinuity set of 
the permeability from a finite set of measurements of the pressure p. Such problems have also been 
studied in the literature. For instance, the paper [46] considers the problem by using multiple level 
set methods in the framework of optimization; and in [27], the authors adopt a Bayesian approach to 
reconstruct the permeability function characterized by layered or channelized structures whose ge¬ 
ometry can be parameterized finite dimensionally. As we consider a finite set of noisy measurements 
of the pressure p, in V*, and the problem is underdetermined and uncertain, the Bayesian approach 
is again natural. We make the significant extension of [27] to consider arbitrary interfaces, requiring 
infinite dimensional parameterization: we introduce a level set parameterization of the domains A> 
as in (2.3) and (2.4). 

Let O : V —> R J denote the collection of J linear functionals on V which are our measurements. 
Because of the estimate (3.5) it is straightforward to see that Q = OoGoFxs bounded as a mapping 
from U into R J and hence that Assumptions 2.1(i) and (iii) hold; it remains to establish (ii). To that 
end, from now on we need slightly higher regularity on /. In particular, we assume that, for some 
q > 2, / G W~ l (L q (D)). Here the space W~\L q (D)) := (W^’ q *(D))* C V* for q* and q 
conjugate: 1/q+l/q* = 1. It is shown in [8] that there exits q 0 > 2 such that the solution of (3.4) 
satisfies 

||Vp||L9(£>) < C\\f\\ w -i( Lq ( D p 

for some C < oo provided 2 < q < q 0 . We assume that such a q is chosen. It then follows that G 
is Lipschitz continuous from L r to V where r := 2q/(q — 2) e [2, oo). To be precise, let p t be the 
solution to the problem (3.4) with Kj, i = 1,2. Then the following is proved in [8]: for any q > 2, 

\\pi — P211< -||Vpi||z,9(u)||/Ci — K 2 IU r (D) 

/^min 


provided Vpi e L q (D). 

Hence G : L r (D) —y V is Lipschitz under our assumption that / e W~ 1 (L q (D)) for some 
q G (2,o). . By viewing F : U —» L r (D), it follows from Proposition (2.4) and Proposition 
(2.5) that Assumptions (2.1) (ii) holds with both Gaussian priors defined in subsection 2.5. As a 
consequence Theorem 2.2 also applies in the groundwater flow model. 
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4 Numerical Experiments 

Application of the theory developed in subsection 2.3 ensures that, for the choices of Gaussian 
priors discussed in subsection 2.5, the posterior measure on the level set is well defined and thus 
suitable for numerical interrogation. In this section we display numerical experiments where we 
characterize the posterior measure by means of sampling with MCMC. In concrete we apply the 
preconditioned Crank-Nicolson (pCN) MCMC method explained in [16]. We start by defining this 
algorithm. Assume that we have a prior Gaussian measure Af(0, C ) on the level set function u and a 
posterior measure n v given by (2.7). Define 

a(u, v ) = min{l, exp(d>(w) — < fi(n))} 

and generate {u^}k >o as follows: 

Algorithm 4.1 (pCN-MCMC). 

Set k = 0 and pick u^> G X 

1. Propose v^ = y/(l — ~^)u^ + /3^ k \ ^ ~ J\f(0,C). 

2. Set = v^ with probability a{u^ k \v^), independently of 

3. Set = u^ otherwise. 

4. k —> k + 1 and return to (i). 

Then the resulting Markov chain is reversible with respect to /j, v and, provided it is ergodic, 
satisfies 

k =0 

for a suitable class of test functions. Furthermore a central limit theorem determines the fluctuations 
around the limit, which are asymptotically of size K~ s. 

4.1 Aim of the experiments 

By means of the MCMC method described above we explore the Bayesian posterior of the level 
set function that we use to parameterize unknown geometry (or discontinuous model parameters) 
in the geometric inverse problems discussed in Section 3. The first experiment of this section con¬ 
cerns the inverse potential problem defined in subsection 3.1. The second and third experiments 
are concerned with the estimation of geologic facies for the groundwater flow model discussed in 
subsection 3.2. The main objective of these experiments is to display the capabilities of the level set 
Bayesian framework to provide an estimate, along with a measure of its uncertainty, of unknown 
discontinuous model parameters in these test models. We recall that for the inverse potential prob¬ 
lem the aim is to estimate the support Di of the indicator function k(x) = (x) that defines the 
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source term of the PDE (3.1) given data/observations from the solution of this PDE. Similarly, given 
data/observations from the solution of the Darcy flow model (3.4), we wish to estimate the inter¬ 
face between geologic facies {A:}” = i corresponding to regions of different structural geology and 
which leads to a discontinuous permeability k(x) = Y^=i K ^o t (%) in the flow model (3.4). In both 
test models, we introduce the level set function merely as an artifact to parameterize the unknown 
geometry (i.e. = {x <E D \ c,_i < u(x) < c*}), or equivalently, the resulting discontinuous field 

k(x). The Bayesian framework applied to this level/ set parameterization then provides us with a 
posterior measure \i y on the level set function u. The push-forward of fj y under the level set map F 
(2.4) results in a distribution on the discontinuous field of interest k. This push-forward of the level 
set posterior F*\i y := \i y o F~ x comprises the statistical solution of the inverse problem which may, 
in turn, be used for practical applications. 

A secondary aim of the experiments is to explore the role of the choice of prior on the posterior. 
Because the prior is placed on the level set function, and not on the model paramerers of direct 
interest, this is a non-trivial question. To be concrete, the posterior depends on the Gaussian prior 
that we put on the level set. While the prior may incorporate our a priori knowledge concerning 
the regularity and the spatial correlation of the unknown geometry (or alternatively, the regions of 
discontinuities in the fields of interest) it is clear that such selection of the prior on the level set may 
have a strong effect on the resulting posterior [i v and the corresponding push-forward F* n y . One of 
the key aspects of the subsequent numerical study is to understand the role of the selection of the 
prior on the level set functions in terms of the quality and efficiency of the solution to the Bayesian 
inverse problem as expressed via the push-forward of the posterior F*/j, y . 

4.2 Implementational aspects 

For the numerical examples of this section we consider synthetic experiments. The PDEs that define 
the forward models of subsection 3 (i.e. expressions (3.1) and (3.4)) are solved numerically, on the 
unit-square, with cell-centered finite differences [4]. In order to avoid inverse crimes [32], for the 
generation of synthetic data we use a much finer grid (size specified below) than the one of size 
80 x 80 used for the inversion via the MCMC method displayed in Algorithm 4.1. 

The Algorithm 4.1 requires, in step (i), sampling of the prior. This is accomplished by param¬ 
eterizing the level set function in terms of the Karhunen-Loeve (KL) expansion associated to the 
prior covariance operator C (See Appendix 2, equation (7.1)). Upon discretization, the number of 
eigenvectors of C equals the dimensions of the discretized physical domain of the model problems 
(i.e. N = 6400 in expression (7.2)). Once the eigendecomposition of C has been conducted, then 
sampling from the prior can be done simply by sampling an i.i.d set of random variables {^.} with 
£i ~ J\f( 0,1) and using it in (7.2). For the present experiments we consider all these KL modes 
without any truncation. However, during the burn-in period (which here is taken to comprise 10 4 it¬ 
erations) of the MCMC method, we find it advantageous to freeze the higher KL modes and conduct 
the sampling only for the lower modes. After the aforementioned burn-in, the sampling is then car- 
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ried out on the full set of KL modes. This mechanism enables the MCMC method to quickly reach 
an “optimal” state where the samples of the level set function provides a field k(x) that is close to 
the truth. However, once this optimal state has been reached, it is essential to conduct the sampling 
on the full spectrum of KL modes to ensure that the MCMC chain mixes properly. More precisely, 
if only the lowest modes are retained for the full chain, the MCMC may collapse into the optimal 
state but without mixing. Thus, while the lowest KL modes determine the main geometric structure 
of the underlying discontinuous field, the highest modes are essential for the proper mixing and thus 
the proper and efficient characterization of the posterior. 

4.3 Inverse Potential Problem 

In this experiment we generate synthetic data by solving (3.1), on a fine grid of size 240 x 240 
with the “true” indicator function /% = I D t displayed in Figure Figure 1 (top). The observation 
operator O = (Oi,..., 0 64 ) is defined in terms of 64 mollified Dirac deltas {Ojjjti centered at 
the measurement locations display as white squares along the boundary of the domain in Figure 1 
(top). Each coordinate of the data is computed by means of (3.3) with p from the solution of the PDE 
with the aforementioned true source term and by adding Gaussian noise rjj with standard deviation 
of 10% of the size of the noise-free measurements (i.e. of Oj (|^)). We reiterate that, in order to 
avoid inverse crimes [32], we use a coarser grid of size 80 x 80 for the inversion via the MCMC 
method (Algorithm 4.1). The parameterization of Lfi in terms of the level set function is given by 
D\ = {x G D | u(x) < 0} (i.e. by simply choosing c 0 = —oo and c\ = 0 in (2.3)). 

For this example we consider a prior covariance C of the form presented in (2.9) for some choices 
of L in the correlation function. We construct C directly from this correlation function and then we 
conduct the eigendecomposition needed for the KL expansion and thus for sampling the prior. In 
Figure 2 we display samples from the prior J\T( 0, C) on the level set function u (first, third and fifth 
rows) and the corresponding indicator function k — I Dl (second, fourth and sixth rows) for (from 
left to right) L = 0.1,0.15,0.2,0.3,0.4. Different values of L in (2.9) clearly result in substantial 
differences in the spatial correlation of the zero level set associated to the samples of the level set 
function. The spatial correlation of the zero level set funtion, under the prior, has significant effect on 
I/jj which we use as the right-hand side (RHS) in problem (3.1) and whose solution, via expression 
(3.3), determines the likelihood (2.6). It then comes as no surprise that the posterior measure on the 
level set is also strong;y dependent on the choice of the prior via the parameter L. We explore this 
effect in the following paragraphs. 

In Figure 3 we present the numerical results from different MCMC chains computed with dif¬ 
ferent priors corresponding to the aforementioned choices of L. The MCMC mean of the level 
set function is displayed in the top row of Figure 3 for the choices (from left to right) L = 
0.1,0.15,0.2,0.3,0.4. We reiterate that although the MCMC method provides the characterization 
of the posterior of the level set function, our primary aim is to identify the field k(x) = Id 1 (x) that 
determines the RHS of (3.1) by means of conditioning the prior Af(0, C ) to noisy data from (3.3). 
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A straightforward estimate of such field can be obtained by mapping, via the level set map (2.4), 
the posterior mean level set function denoted by u into the corresponding field F(u(x )) = k(x) = 
I-/j ] ( x ) where D\ = {x e D \u(x) < 0}. We display k (x) = Ijj ( x ) in the top-middle row of Figure 
3 along with the plot of the true field F = I D t (right column) for comparison. 

As mentioned earlier, we are additionally interested in the push-forward of the posterior measure 
of the level set function u under the level set map (i.e. ( F*/x y )(du )). We characterize F*/i y by 
mapping under F our MCMC samples from /i y . In Figure 3 we present the mean (bottom-middle) 
and the variance (bottom) of F*/j, y . Figure 4 shows some posterior (MCMC) samples u of the 
level set function (first, third and fifth rows) and the corresponding level set map F(u) = I d 1 with 
D\ = {x G D | u(x) < 0} associated to these posterior samples (second, fourth and sixth rows). 

The push-forward of the posterior measure under the level set map (i.e. F*/x y ) thus provides 
a probabilistic description of the inverse problem of identifying the true F = I D t. We can see 
from Figure 3 that, for some choices of L, the mean of F* /j y provides reasonable estimates of the 
truth. However, the main advantage of the Bayesian approach proposed here is that a measure of 
the uncertainty of such estimate is also obtained from F*n y . The variance (Figure 3 bottom), for 
example, is a measure of the uncertainty in the location of the interface between the two regions D 
and D \ D\. 

The results displayed in Figure 3 show the strong effect that the selection of the prior has on the 
posterior measure /A and the corresponding pushforward measure F*fi y . In particular, there seems 
to be a critical value L = 0.2 above of which the corresponding posterior mean on F*fj y provides 
a reasonable identification of the true II D t with relatively small variance. This critical value seems 
to be related to the size and the shape of the inclusions that determines the true region D\ (Figure 
1 (top)). It is intuitive that posterior samples that result from very small spatial correlation cannot 
easily characterize these inclusions accurately unless the data is overwhelmingly informative. The 
lack of a proper characterization of the geometry from priors associated with small L is also reflected 
with larger variance around the interface. It is then clear that the capability of the proposed level 
set Bayesian framework to properly identify a shape D\ (or alternatively its indicator function I D t) 
depends on properly incorporating, via the prior measure, a priori information on the regularity and 
spatial correlation of the unknown geometry of D\. 

Since the selection of the prior has such a clear effect on the posterior, it comes as no surprise 
that it also affects the efficiency of the MCMC method as we now discuss. In the bottom-right panel 
of Figure 1 we show the autocorrelation function (ACF) of the first KL mode of the level set func¬ 
tion from different MCMC chains with different priors corresponding to our different choices of 
correlation length L in (2.9). The tunable parameters in the pCN-MCMC method are fixed for these 
experiments. We recall from Figure 3 that larger values of L result in a mean level set whose cor¬ 
responding indicator function better captures the spatial structures form the truth and with smaller 
variance around the interface. However, the larger the value of L the slower the decay of the ACF. 
From these ACF plots, we note that even for the apparent optimal value of L = 0.3, our MCMC 
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method produces samples that are highly correlated and thus very long chains may be needed in 
order to produce a reasonable number of uncorrelated samples needed for statistical analysis. For 
this particular choice of L — 0.3 we have conducted 50 multiple MCMC chains of length 10 6 (after 
burn-in period) initialized from random samples from the prior. In Figure 1 (bottom-left) we show 
the Potential Scale Reduction Factor (PSRF) [10] computed from MCMC samples of the level set 
function (red-solid line) and the corresponding samples under F (i.e. the Ip/s) (blue-dotted line). 
We observe that the PSRF goes below 1.1 after (often taken as an approximate indication of conver¬ 
gence [10]); thus the Gelman-Rubin diagnostic [10] based on the PSRF is passed for this selection 
of L. The generation of multiple independent MCMC chains that are statistically consistent opens 
up the possibility of using high-performance computing to enhance our capabilities of properly ex¬ 
ploring the posterior. While we use a relatively small number of chains as a proof-of-concept, the 
MCMC chains are fully independent and so the computational cost of running multiple chains scales 
with the number of available processors. 

The 5 x 10 7 samples that we obtained from the 50 MCMC chains are combined to provide a 
full characterization of the posterior /F on the level set and the corresponding push-forward F*fF 
(i.e. The I^’s computed from D\ with posterior samples u). We reemphasize that our aim is the 
statistical identification of . Therefore, in order to obtain a quantity from the true I D t against 
to which compare the computed push-forward of the level set posterior, we consider the Discrete 
Cosine Transform (DCT) of the true field I D . Other representations/expansions of the true field 
could be considered for the sake of assessing the uncertainty of our estimates with respect to the 
truth. In Figure 5 we show the prior and posterior densities of the first DCT coefficients of 
where Di = {x e D \u(x) < 0} with u from our MCMC samples (the vertical dotted line 
corresponds to the DCT coefficient of the true Ipt). We can observe how the push forward posterior 
are concentrated around the true values. It is then clear how the data provide a strong conditioning 
on the first DCT coefficients of the discontinuous field that we aim at obtaining with our Bayesian 
level set approach. 
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Figure 1: Inverse Potential. Top: True source term k) = I £) \ of eq. (3.1). Bottom-left: PSRF from 
multiple chains with L = 0.3 in (2.9). Bottom-right: ACF of first KL mode of the level set function 
from single-chain MCMC with different choices of L. 
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Figure 2: Inverse Potential. Samples from the prior on the level set function u (first, third and fifth 
rows) for (from left to right) L = 0.1,0.15,0.2,0.3,0.4. Corresponding I Dl with D\ = {x G 
D | u(x) < 0} (second, fourth and sixth rows) associated to each of these samples from the level set 
function. 
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Figure 3: Inverse Potential. MCMC results for (from left to right) L = 0.1,0.15, 0.2, 0.3,0.4 in the 
eq. (2.9). Top: Posterior mean level set function u (computed via MCMC). Top-middle: Plot of % | 
with Di = {x e D \u(x) < 0} (the truth I D t is presented in the right column). Bottom-middle: 
Mean of lr> 1 where D\ = {x G D \u(x) < 0} and u ’s are the posterior MCMC samples (the truth 
is presented in the right column). Bottom: Variance of lo, where D\ = (:r £ D \u(x) < 0} and 
u’s are the posterior MCMC samples 





















Figure 4: Inverse Potential. Samples from the posterior on the level set u (first, third and fifth rows) 
for (from left to right) L = 0.1,0.15, 0.2, 0.3,0.4. Corresponding l Dl where Di = {x e D \u(x) < 
0} (second, fourth and sixth rows) associated to each of these samples from the level set function. 
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Figure 5: Inverse Potential. Densities of the prior and posterior of various DCT coefficients of the 
Ii) 1 where D y = {x e D |u(x) < 0} obtained from MCMC samples on the level set u for L = 0.3 
(vertical dotted line indicates the truth). 
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4.4 Structural Geology: Channel Model 

In this section we consider the inverse problem discussed in subsection 3.2. We consider the Darcy 
model (3.4) but with a more realistic set of boundary conditions that consist of a mixed Neumman 
and Dirichlet conditions. For the concrete set of boundary conditions as well as the right-hand-side 
we use for the present example we refer the reader to [29, Section 4]. This flow model, initially used 
in the seminal paper of [14], has been used as a benchmark for inverse problems in [29, 22, 24], 
While the mathematical analysis of subsection is 3.2 conducted on a model with Dirichlet boundary 
conditions, in order to streamline the presentation, the corresponding extension the case of mixed 
boundary conditions can be carried out with similar techniques. 

We recall that the aim is to estimate the interface between regions Di of different structural geol¬ 
ogy which result in a discontinuous permeability n of the form (2.2). In order to generate synthetic 
data, we consider a true k\x) = ^? =1 K, 7 II D t with prescribed (known) values of k i = 7, k 2 = 50 
and k 3 = 500. This permeability, whose plot is displayed in Figure 6 (top), is used in (3.4) to gen¬ 
erate synthetic data collected from interior measurement locations (white squares in Figure 6). The 
estimation of k is conducted given observations of the solution of the Darcy model (3.4). To be con¬ 
crete, the observation operator O = (Oi,..., 0 2 5 ) is defined in terms of 25 mollified Dirac deltas 
centered at the aforementioned measurement locations and acting on the solution p of the 
Darcy flow model. For the generation of synthetic data we use a grid of 160 x 160 which, in order 
to avoid inverse crimes [32], is finer than the one used for the inversion (80 x 80). As before, ob¬ 
servations are corrupted with Gaussian noise proportional to the size of the noise-free observations 
( Oj(p ) in this case). 

For the estimation of k with the proposed Bayesian framework we assume that knowledge of 
three nested regions is available with the permeability values {Kj}f = 1 that we use to define the 
true k 1 . Again, we are interested in the realistic case where the rock types of the formation are 
known from geologic data but the location of the interface between these rocks is uncertain. In 
other words, the unknowns are the geologic facies D t that we parameterize in terms of a level set 
function, i.e. Di = {x G D |q_i < u(x) < q} with c 0 = — 00 , C\ = 0, c 2 = 1, c 3 = 00 . 
Similar to the previous example, we use a prior of the form (2.9) for the level set function. In Figure 
7 we display samples from the prior on the level set function (first, third and fifth rows) and the 
corresponding permeability mapping under the level set map (2.4) F{u){x) = k(x') = Y^t=\ k >^d, 
(second, fourth and sixth rows) for (from left to right) L = 0.2, 0.3, 0.35, 0.4,0.5. As before, we 
note that the spatial correlation of the covariance function has a significant effect on the spatial 
correlation of the interface between the regions that define the interface between the geologic facies 
(alternatively, the discontinuities of k). Longer values of L provide k’s that seem more visually 
consistent with the truth. The results from Figure 8 show MCMC results from experiments with 
different priors corresponding to the aforementioned choices of L. The posterior mean level set 
function u is displayed in the top row of Figure 8. The corresponding mapping under the level set 
function k = J2t=i K ^o, (with D, = (a; <S D |q_i < u(x) < c,}) is shown in the top-middle. 
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Similar to our discussion of the preceding subsection, for the present example we are interested 
in the push-forward of the posterior [i y under the level set map F. More precisely, F*n y provides 
a probability description of the solution to the inverse problem of finding the permeability given 
observations from the Darcy flow model. In Figure 8 we present the mean (bottom-middle) and the 
variance (bottom) of F(n y ) characterized by posterior samples on the level set function mapped 
under F. In other words, these are the mean and variance from the k’s obtained from the MCMC 
samples of the level / set function. As in the previous example, there is a critical value of L = 
0.3 below of which the posterior estimates cannot accurately identify the main spatial features of 

Figure 9 shows posterior samples of the level set function (first, third and fifth rows) and the 
corresponding k (second, fourth and sixth rows). The posterior samples, for values of L above the 
critical value L = 0.3, capture the main spatial features from the truth. There is, however, substantial 
uncertainty in the location of the interfaces. Our results offer evidence that this uncertainty can 
be properly captured with our level set Bayesian framework. Statistical measures of F*/j, y (i.e. 
the posterior permeability measure on k) is essential in practice. The proper quantification of the 
uncertainty in the unknown geologic facies is vital for the proper assessment of the environmental 
impact in applications such as C0 2 capture and storage, nuclear waste disposal and enhanced oil 
recovery. 

In Figure 6 (bottom-right) we show the ACF of the first KL mode of level set function from 
different MCMC chains corresponding to different priors defined by the choices of L indicated pre¬ 
viously. In contrast to the previous example, here we cannot appreciate substantial differences in the 
efficiency of the chain with respect to the selected values of L. However, we note that ACF exhibits 
a slow decay and thus long chains and/or multiple chains are need to properly explore the poste¬ 
rior. For the choice of L = 0.4 we consider 50 multiple MCMC chains. Our MCMC chains pass 
the Gelman-Rubin test [10] as we can note from Figure 6 (bottom-left) where we show the PSRF 
computed from MCMC samples of the level set function u (red-solid line) and the corresponding 
mapping, under the level set map, into the permeabilities k (blue-dotted line). As indicated earlier, 
we may potentially increase the number of multiple chains and thus the number of uncorrelated 
samples form the posterior. 

Finally, in Figure 10 we show the prior and posterior densities of the first DCT coefficients on the 
k obtained from the MCMC samples of the level set function (the vertical dotted line corresponds 
to the DCT coefficient of the truth /-d). For some of these modes we clearly see that the posterior is 
concentrated around the truth. However, for the mode Cm we note that the posterior is quite close 
to the prior indicating that the data have not informed this mode in any significant way. 

4.5 Structural Geology: Layer Model 

In this experiment we consider the groundwater model (3.4) with the same domain and measurement 
configurations from the preceding subsection. However, for this case we define the true permeability 

displayed in Figure 11 (top). The permeability values are as before. The generation of synthetic 
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Figure 6: Identification of structural geology (channel model). Top: True k in eq. (3.4). Bottom-left: 
PSRF from multiple chains with L = 0.4 in (2.8). Bottom-right: ACF of first KL mode of the level 
set function from single-chain MCMC with different choices of L. 


data is conducted as described in the preceding subsection. For this example we consider the Gaus¬ 
sian prior on the level set defined by (2.8). Since for this case the operator C is diagonalisable by 
cosine functions, we use the fast Fourier transform to sample from the corresponding Gaussian 
measure Af( 0, C) required by the pCN-MCMC algorithm. 

The tunable parameter a in the covariance operator (2.8) determines the regularity of the cor¬ 
responding samples of the Gaussian prior (see for example [45]). Indeed, in Figure 12 we show 
samples from the prior on the level set function (first, third and fifth rows) and the corresponding k 
(second, fourth and sixth rows) for (from left to right) a = 1.5, 2.0, 2.5, 3.0, 3.5. Indeed, changes in 
a have a dramatic effect on the regularity of the interface between the different regions. We therefore 
expect strong effect on the resulting posterior on the level set and thus on the permeability. 

In Figure 13 we display numerical results from MCMC chains with different priors correspond¬ 
ing to (from left to right) a = 1.5,2.0, 2.5,3.0,3.5. In Figure 13 we present the MCMC mean of the 
level set function.The corresponding k is shown in the top-middle of Figure 13. In this figure we 
additionally display the mean (bottom-middle) and the variance (bottom) of the Ac’s obtained from 
the MCMC samples of the level set function. Above a critical value a = 2.5 we obtain a reasonable 
identification of the layer permeability with a small uncertainty (quantified by the variance). Figure 
14 shows posterior (MCMC) samples of the level set function (first, third and fifth rows) and the 
corresponding k (second, fourth and sixth rows) for the aforementioned choices of a. 


























A Bayesian Level Set Method 


23 




» 


• » 


i umu 




Figure 7: Identification of structural geology (channel model). Samples from the prior on the level 
set (first, third and fifth rows) for (from left to right) L = 0.2,0.3,0.35,0.4,0.5. Pushforward onto 
k (second, fourth and sixth rows) associated to each of these samples from the level set function. 



































I I I 


24 


M. Iglesias, Y. Lu and A. Stuart 



Figure 8: Identification of structural geology (channel model). MCMC results for (from left to right) 
L = 0.2,0.3,0.35,0.4,0.5 in the eq. (2.8). Top: MCMC mean of the level set function. Top-middle: 
k associated to the mean of the level set function (true k is displayed in the last column). Bottom- 
middle: Mean of the n. Bottom: Variance of k 
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Figure 9: Identification of structural geology (channel model). Samples from the posterior on the 
level set (first, third and fifth rows) for (from left to right) L = 0.2,0.15,0.2,0.3,0.4. log(/c) (second, 
fourth and sixth rows) associated to each of these samples from the level set function. 
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Figure 10: Identification of structural geology (channel model). Densities of the prior and posterior 
of some DCT coefficients of the Ac’s obtained from MCMC samples on the level set for L — 0.4 
(vertical dotted line indicates the truth). 


Figure 11 (bottom-right) shows the ACF of the first KL mode of level set function from MCMC 
experiments with different priors with cc’s as before. The efficiency of the MCMC chain does not 
seem to vary significantly for the values above the critical value of a. However, as in the previous 
examples a slow decay in the ACF is obtained. An experiment using 50 multiple MCMC chains 
initialized randomly from the prior reveals that the Gelman-Rubin diagnostic test [10] is passed for 
a = 2.5 as we can observe from Figure 11 (bottom-left) where we the display PSRF from MCMC 
samples of the level set function (red-solid line) and the corresponding mapping into the k (blue- 
dotted line). In Figure 15 we show the prior and posterior densities of the DCT coefficients on the k 
obtained from the MCMC samples of the level set function (the vertical dotted line corresponds to 
the truth DCT coefficient). We see clearly that the DCT coefficients are substantially informed by 
the data although the spread around the truth confirms the variability in the location of the interface 
between the layers that we can ascertain from the posterior samples (see Figure 14). 

5 Conclusions 

The primary contributions of this paper are: 

• We have formulated geometric inverse problems within the Bayesian framework. 

• This framework leads to a well-posedness of the level set approach, something that is hard to 
obtain in the context of classical regularization techniques for level set inversion of interfaces. 
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Figure 11: Identification of structural geology (layer model). Top: True k in eq. (3.4). Bottom-left: 
PSRF from multiple chains with a = 2.5 in (2.9). Bottom-right: ACF of first KL mode of the level 
set function from single-chain single-chain MCMC with different choices of a. 


• The framework also leads to the use of state-of-the-art function-space MCMC methods for 
sampling of the posterior distribution on the level set function. An explicit motion law for the 
interfaces is not needed: the MCMC accept-reject mechanism implicitly moves them. 

• We have studied two examples: inverse source reconstruction and an inverse conductivity 
problem. In both cases we have demonstrated that, with appropriate choice of priors, the 
abstract theory applies. We have also highlighted the behavior of the algorithms when applied 
to these problems. 

• The fact that no explicit level set equation is required helps to reduce potential issues arising 
in level set inversion, such as reinitialization. In most computational level set approaches [13], 
the motion of the interface by means of the standard level set equation often produces level 
set functions that are quite flat. For the mean curvature flow problem, such fattening phenom¬ 
ena is firstly observed in [21] where the surface evolution starts from a ’’figure eight” shaped 
initial surface; in addition it has been shown to happen even if the initial surface is smooth 
[42], This causes stagnation as the interface then move slowly. Ad-hoc approaches, such as 
redistancing/reinitializing the level set function with a signed distance function, are then em¬ 
ployed to restore the motion of the interface. In the proposed computational framework, not 
only does the MCMC accept-reject mechanism induce the motion of the intertace, but it does 
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Figure 12: Identification of structural geology (layer model). Samples from the prior on the level 
set (first, third and fifth rows) for (from left to right) a = 1.5,2.0,2.5,3.0,3.5 in (2.9). k (second, 
fourth and sixth rows) associated to each of these samples from the level set function. 
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Figure 13: Identification of structural geology (layer model). MCMC results for (from left to right) 
a = 1.5,2.0,2.5,3.0,3.5 in the eq. (2.9). Top: MCMC mean of the level set function. Top-middle: 
k associated to the mean of the level set function (true k is displayed in the last column). Bottom- 
middle: Mean of the n. Bottom: Variance of k 
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Figure 14: Identification of structural geology (layer model). Samples from the posterior on the level 
set (first, third and fifth rows) for (from left to right) a = 1.5,2.0,2.5,3.0,3.5 in the eq. (2.9). k 
(second, fourth and sixth rows) associated to each of these samples from the level set function. 
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Figure 15: Identification of structural geology (layer model). Densities of the prior and posterior 
of some DCT coefficients of the Ac’s obtained from MCMC samples on the level set for L — 0.4 
(vertical dotted line indicates the truth). 


so in a way that avoids creating flat level set functions underlying the permeability. Indeed, 
we note that the posterior samples of the level set function inherit the same properties from 
the ones of the prior. Therefore, the probability of obtaining a level set function which takes 
any given value on a set of positive measure is zero under the posterior, as it is under the prior. 
This fact promotes very desirable, and automatic, algorithmic robustness. 

Natural directions for future research include the following: 

• The numerical results for the two examples that we consider demonstrate the sensitive de¬ 
pendence of the posterior distribution on the length-scale parameter of our Gaussian priors. It 
would be natural to study automatic selection techniques for this parameter, including hierar¬ 
chical Bayesian modelling. 

• We have assumed that the values of ACj on each unknown domain D t are both known and 
constant. It would be interesting, and possible, to relax either or both of these assumptions, as 
was done in the finite geometric parameterizations considered in [27]. 

• The numerical results also indicate that initialization of the MCMC method for the level set 
function can have significant impact on the performance of the inversion technique; it would 
be interesting to study this issue more systematically. 

• The level set formulation we use here, with a single level set function and possibly multiple 
level set values c t has been used for modeling island dynamics [38] where a nested structure 
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is assumed for the regions {A}"=i see Figure 16(a). However, we comment that there exist 
objects with non-nested regions, such as the one depicted in Figure 16(b), which can not be 
represented by a single level set function. . It would be of interest to extend this work to the 
consideration of vector-valued level set functions. In the case of binary obstacles, it is enough 
to represent the shape via a single level set function (cf. [43]). However, in the case of n-ary 
obstacles or even more complex geometric objects, the representation is more complicated; 
see the review papers [18, 19, 46] for more details. 



Figure 16: Nested regions and non-nested regions 


• The Bayesian framework could be potentially combined with other parameterizations of un- 
kwon geometries. For example, the pluri Gaussian approach has been used with EnKF in [35] 
to identify geologic facies. 
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6 Appendix 1 

Proof of Theorem 2.2. Notice that the random variable y\u is distributed according to the measure 
Q u , which is the translate of Q 0 by Q(u), satisfying Q u -C Qo with Radon-Nikodym derivative 


dQo 


(y) oc exp(-$(w;?/)); 


where $ : U x R J —>■ E is the least squares function given in (2.6). Thus for the given data y, 
<b(u\ y) is up to a constant, the negative log likelihood. We denote by u 0 the product measure 


u 0 (du,dy ) = p 0 (du)Q 0 (dy). 


( 6 . 1 ) 


Suppose that $(•,•) is u 0 measurable, then the random variable (u,y) G U x Y is distributed 
according to u(du, dy) with 

u,y) OC exp(-$(u;y)). 

From Assumptions 2.1(ii) and the continuity of $(m; y) with respect to y, we know that $(•; •) : 
U x Y — » M is continuous ^-almost surely. Then it follows from Lemma 6.1 below that $(•; •) 
is z'o-measurable. On the other hand, by Assumptions 2.1(i), for | vy | p < r, we obtain the upper and 
lower bound for Z, 

0 < exp(—A(r, « min )) <Z= exp (-^(u;y))po(du) < 1 

Ju 

Thus the measure is normalizable and applying the Bayes’ Theorem 3.4 from [17] yields the exis¬ 
tence of p y . 

Let Z = Z{y) and Z' = Z(y') be the normalization constants for p y and g 11 ' , i.e. 

Z = \ ex-p(-$(u-y))g, 0 (du), Z'= exp (-§(ir,y'))g 0 (du) 


Ju 

We have seen above that 


<u 


1 > Z, Z' > exp (~K(r, K min )) > 0. 

From Assumptions 2. l(iii), 

\Z-Z'\< J \exp(-$>(u]y))-exp(-$(u]y'))\g 0 (du) < J \$(u]y)-$(u]y')\g 0 (du) < C(r)\y-y' 

Thus, by the definition of Hellinger distance, we have 

2d nen (g y ,g y ') 2 = J ^Z _1/2 exp - (Z') _1/2 exp j/)^ g 0 (du) < h+I 2 
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where 

j - exp ^-^$(u-y' 

h = 2|Z“ 1/2 - (Z')- 1/2 \ 2 J exp(-$(u;y'))ia 0 (du) 

Applying (i) and (iii) in Assumptions 2.1 again yields 

< C(r)\y-y'\ 2 v 

and 

h < C(r)\Z~ 1/2 - (Z')“ V T < C(r)(Z~ 3 V ( Z')~ 3 )\Z - Z'\ 2 < C{r)\y - y'\ 2 v 

Therefore the proof that the measure is Lipschitz is completed by combining the preceding esti¬ 
mates. The final statement follows as in [45], after noting that / G L 2 o ((7; S) implies that / G 
Lly{U-S), since $(•;?/) > 0. □ 

Lemma 6.1. Let U be a separable Banach space and ( U , E, //) be a complete probability space with 
(j-algebra E. If a functional T :£/—>■ M is continuous //-almost surely, i.e. = 1 where M 

denotes the set of the continuity points of J then T is E-measurable. 

Proof. By the definition of measurability, it suffices to show that for any c > 0 

S := {u G U | F{u) > c} G E. 

One can write S as S = (5'nA7)U(S'\ M). Since T is continuous //-almost surely, M is measurable 
and n(M) = 1. It follows from the completeness of the measure // that S\M is measurable and 
y(S \ M ) = 0. Now we claim that S fl M is also measurable. Denote Bs(u) C U to be the ball of 
radius 8 centered at u G U. For each v G S fl M, as J 7 is continuous at v, there exists 8 V > 0 such 
that if v' G Bs v (v), then J-(v') > c. Therefore S H M can be written as 

SDM = Mf] U B Su (//) 

vesnM 

that is the intersection of the measurable set M with the open set U^esriA/ B Sv (y ). So S n M is 
measurable. Then it follows that T is E-measurable. □ 

Proof of Proposition 2.4. Let {// s } denote any approximating family of level set functions 

with limit u as e — > 0 in C(D; M) : \\u e — u\\ c ^ < e —> 0. Let D i e , D { [ £ be the sets defined in (2.3) 
associated with the approximating level set function u £ and define k = F(u) by (2.4) and, similarly, 
n e := F{u e ). Let rn(A) denote the Lebesgue measure of the set A. 

Suppose that m(Df) = 0, i = 1, • • • ,n — 1. Let {ip} be the above approximating functions. We 
shall prove \\k £ — k\\l<i(d) —>■ 0. In fact, we can write 

«e0*0 - «(*) = E"=i E"=i(«i - ^D^nDjix) 

= ~ K j)^D i>£ nDj(x). 
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By the definition of u e , for any x £ D 

u(x) — £ < U e {x) < u(x) + £ ( 6 . 2 ) 

Thus for | j — i\ > 1 and £ sufficiently small, D t e n Dj = 0. For the case that \i — j\ = 1, from 
( 6 . 2 ), it is easy to verify that 


A,£ n D i+ 1 C D i e := {x £ D I a < u(x) < Ci + £} ->• £>°> * = , n - 1 (6.3) 

A,e n A -1 c A-i,e := {xefi Q-i - £ < «(x) < c;_i} ->0, i = 2, • • • n (6.4) 


as £ —» 0. By this and the assumption that m{Df) = 0, we have that rn(D l £ D Dj) —> 0 if i ^ j. 
Furthermore, the Lebesgue’s dominated convergence theorem yields 

n p 

||k £ tt|| Li(D) ^ y / I K'i Kjl^diE > 0 

• • 1 • / • A,- e n£>, 

as £ —> 0. Therefore, F is continuous at u. 

“=>” We prove this by contradiction. Suppose that there exists i* such that m(D%) ± 0. We 
define u £ := u — e, then it is clear that ||tt e — u\\ c ^ —>■ 0 as £ —>■ 0. By the same argument used in 
proving the sufficiency, 


\k f — k 


\ q 

\L*(D) 


n— 1 r 

£/. - 

^—1 J Di j£ UDi £ 


\K i+ i~Ki\ q dx ->• 


n— 1 « 

J] / l Ki + r 


-Ki\ q dx > / |Kj* + i—| 9 dx > 0 

Jd° 


where we have used m(Df») 7 ^ 0 in the last inequality. However, this contradicts with the assump¬ 


tion that F is continuous at u. 


□ 


7 Appendix 2 

Recall the Gaussian measure /j 0 = 7V"(0. C) on the function space F where C = Ci,H = Hi,i = 1,2 
given in subsection 2.5. These measures can be constructed as Gaussian random fields. 

Let (f2, JF, P) be a complete probability space, i.e. if A £ JF with P(^4) = 0, then P(F) = 0 
for B C A. A random field on D is a measurable mapping u : D x fl — » M. Thus, for any 
x £ D, u{x\ •) is a random variable in M; whilst for any ui £ (2, «(•; w) : D —» M is a vector field. 
Denote by (M N , £>(M N ), P) the probability space of i.i.d standard Gaussian sequences equipped with 
product cr-algebra and product measure. In this paper, we consider (O, JF, P) as the completion 
of (M N , £>(M N ), P) in which case the cr-algebra JF consists of all sets of the type Au B, where 
A £ F(M N ) and B C N £ F(M N ) with P(iV) = 0. Let uj = {£&} Ai £ ^ = IR N be an i.i.d sequence 
with ~ A/"(0,1), and consider the random function u £ H defined via the Karhunen-Loeve 


expansion 
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u(x\u) = T(u) := Y, y/\k£k<f>k(x), (7.1) 

k =1 

where (A/,,. <A/ c }/A=i is the eigensystem of C. By the theory of Karhunen-Loeve expansions [7], the 
law of the random function u is identical to // () . Recalling that a > 1 , the eigenvalues {A/.} of C\ 
decay like k~ a in two dimensions; whilst the eigenvalues of C 2 will decay exponentially. Moreover, 
we assume further that <pi~ G U and that sup fc Halloo < oo which holds in simple geometries. Due 
to the decaying properties of the eigenvalues of C, the truncated sum 

N 

Tn(u) = ^ \Z^k^k4>k (7.2) 

k =1 

admits a limit T in I/jj>(12; H). By the Kolmogorov Continuity Theorem [17], T is in fact Holder 
continuous P—almost surely; in particular, T G U P-almost surely. Then by Theorem 3.1.2 in [ ], 
we have 7jv —> T in the uniform norm of U, P-almost surely. Since for any N G , 7jv : (12, <&) —>■ 
(i U , 3§(U)) is continuous and thus measurable, we know from the completeness of (12, ,A, P) that 
the limit T is also measurable from ( 12 , &) to (U , &(U)) (see p30 in [44]). The measurability of T 
enables us to define a new measure on (U,&(U)) which we still denote by // 0 by the following: 

/io(A) =P(r _1 (^)) = P({cc G 12 I u(-;cc)gA}) fotAe&(U). (7.3) 

Thus /i 0 is indeed the push-forward measure of P through T. By definition, it is not hard to verify 
that fj, 0 is the Gaussian measure Af(0, C ) on (U. 3§(U)). In addition, suppose that B c N <E &(U) 
with Ho(N) = 0; if we still define // 0 ( B ) according to (7.3), then // 0 ( B ) = P (T~ 1 (B)) = 0 by 
the fact that T~ 1 (B) c T _l (N) and the completeness of (12, , P). Denote by E the smallest a 

algebra containing 3S(U) and all sets of zero measure under // 0 so that any set E G E is of the form 
E = A U B, where A G &{U) and B <Z N E &(U) with /j, 0 (N) = 0. Then ( U , E, /j, 0 ) is complete. 

We comment that although a Gaussian measure is usually defined as a Borel measure in the liter¬ 
ature (see e.g. [7]), it is more convenient to work with a complete Gaussian measure in this paper; in 
particular, the completeness of /i 0 is employed to show the measurability of the observational map 
in level set based inverse problems. 

Considering a Gaussian random function u(-;u) with co G 12, for any level constant c G M, we 
define the random level set 

D° c = D° c (u{-]uj)) = D° c {<jj) := {x | u(x\uj) = c}. (7.4) 

Recall that the measure space (U, E, /j 0 ) is the push-forward of (12, P) under T. We define the 
functional M c : U —> M by 

A i c u = m({x | u(x) = c}) 
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and the composition 1Z C := M. c o T, as illustrated in the following commutative diagram: 

(12,#", p) — > (u,n,no) 



(R,@(R)) 


Lemma 7.1. For any c E R, M c is E-measurable and 7 Z c is ^-measurable so that m(D° c ) is a 
random variable on both (U, E, /i 0 ) and (12, P). 

Proof. To prove A4 C is E-measurable, we only need to check the set A t : = (u E U M c u > t} E E 
for any t 6 1. Since Af c is a non-negative map, for t < 0, it is obvious that A t = U and hence 
measurable. Now we claim that A t is closed in U for t > 0. To that end, let {u n }fL x be a sequence 
of functions in A t such that | u n — u \ \ u —>■ 0 for some u E U as n —>■ oo. We prove that u E A t . Since 
\\u n ~ u \\u —■> 0, there exists a subsequence which is still denoted by u n such that | u n — u \\u < 1 /n. 
By the definition of A t , u n e A t means that m({x E D \ u n (x) = c}) > t for all n. Moreover, from 
the construction of u n , {x E D \ u n (x) = c} C B n := {x E D \ \ u(x) — c\ < 1/n}, which implies 
that m(B n ) > t. Noting that 

{x E D | u(x) = c] = fl ™ =1 B n 

and that B n is decreasing, we can conclude that m({x E D | u(x) = 0}) > t, i.e. u E A t . So A t is 
closed for t > 0. Then it follows from the measurability of T that 7 Z c is ^"-measurable. Therefore 
m{D ° c ) is a random variable on both (U, E, fi Q ) and (12, &, P). □ 

The following theorem demonstrates that m{D ° c ) vanishes almost surely on both measure spaces 
above. 

Proposition 7.2. Consider a random function u drawn from one of the Gaussian probability mea¬ 
sures go on (27, E) given in subsection 2.5. For c 6 l, the random level set D° c of u is defined by 
(7.4). Then 

(i) m{D° c ) = 0,P-almost surely; 

(ii) m(D °) = 0 , /iq -almost surely. 

Proof, (i) For any fixed x E D, since the point evaluation u(x) acts as a bounded linear functional 
on U, u(x] •) is a real valued Gaussian random variable, which implies P({w | u(x,cu) = c }) = 0. 
Moreover, noting that the random field u : D x Cl —>Misa measurable map, if we view rn(Df) as 
a random variable on 12 , then 


E[m(D° c )}= / m (D°(u)) dP(u) = / / !{* | u (x-,w)=c} dx dP(w) 



Fubini 


a JR d 


I{(a:,a>)|«(a:;a>)=c} dxdP(a;) = 



R d JQ 


| u(x;ui)=c} dP(tu) dx 


Cl JR d 

^{(x, w) | u(x;ui)=c} dP((j) dx 

/R d JCI 

P({cu | u(x; oj) = c}) dx = 0 
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Noting that m(D ° c ) > 0, we obtain m{D° c ) = 0, P-almost surely. 

(ii) Recall that A t = {u G U \ A4 c u > t} defined in Lemma 7.1 is closed in U for any t > 0. 
Thus the set A := {u e U | m({x | u(x) = c}) = 0} = ( Uj? =1 Ai/ k ) c = n ^ =1 A\, k is a Borel set of 
U and measurable. Since /i 0 is the push-forward measure of P under T, 

= P (T-\A)) = P({w | m(B»(w)) = 0}) = 1 

where the last equality follows from (i). □ 


